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1 Introduction 



The XXZ Heisenberg spin i finite chain of length M ||, in a constant magnetic field 
parallel to the anisotropy direction z, is defined by the Hamiltonian, 

HxXZ = {o-m^m+l + + A ( CT m cr m+l ~ 1) ~ 2°™}' 

m=l 

where a m , a = x,y,z, are the Pauli matrices associated to the site m of the chain and acting 
in the two-dimensional space TL m at site m, A is the anisotropy parameter and h is a constant 
magnetic field in the z direction. 

The eigenstates and the spectrum of this model have been described using the Bethe ansatz 
method in [§,0], and later by its algebraic version in ||. 

The aim of this article is to describe the n-point correlation functions of the local spin 
operators of this model in the thermodynamic limit, in the presence of a constant magnetic field 
h, for the anisotropy parameter A > — 1 and at zero temperature. If the anisotropy parameter is 
less than —1, the ground state of the Hamiltonian ( |1.1|) is ferromagnetic and all such correlation 
functions can be calculated easily. In fact let us denote by | ip g } the ground state in the massless 
regime and any one of two ground states constructed by the algebraic Bethe ansatz in the 

massive regime. Let E m n ' <Lm be the elementary operators acting on 7i m at site m as the 2x2 

, m £ / e , 

matrices = 5[ e >5k e , and Y\ Ej' J an Y product of such elementary operators from site one 

3=1 

to to. We will compute the following correlation functions, 

(<Mn4'''T* s } 

F - ({ ^»- »M ■ (1 ' 2) 

Any n-point correlation function can be obtained from these building blocks. 

Our approach is based on the algebraic Bethe ansatz formulation of this model || and 
the actual resolution of its associated quantum inverse scattering problem for local spins given 
in ||. The first results in the thermodynamic limit were obtained using this method for the 
spontaneous magnetization of the XXZ model in Q. 

Previous results about correlation functions were obtained in || and || . In the case h = our 
results agree with the one following the q-deformed Knizhnik-Zamolodchikov (KZ) equations Q 
in both the massless (—1 < A < 1) [^] and massive (A > 1) [|TT) regimes. 

Before describing our method, let us first briefly recall the problem to be solved, the main 
difficulties to overcome, and the approaches that have been used previously in this context. 

In fact, the computation of exact correlation functions is one of the most challenging and 
longstanding problem in the domain of integrable models of field theory and statistical mechan- 
ics. This is a central question both from the fundamental and mathematical point of view, and 
also when one is willing to apply such models for example to condensed matter phenomena. 

Despite the great amount of remarkable works (see |8, 12, 13 1 and references therein), no 



general method for computing exact and explicit manageable expressions for the correlation 
functions of generic integrable models is available at present. 

In this context, the XXZ Heisenberg spin-| chain in the anti-ferromagnetic regime has already 
been used as a representative example of integrable model, with non-trivial interaction and 
bound states (in particular it is not equivalent to free fermions, except at A = 0), associated to 
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a q-deformed trigonometric ii-matrix, and which contains all the basic difficulties in computing 
correlation functions. 

In fact, while the ground state of such models in the thermodynamic limit can be described 
very precisely, for example by means of the Bethe ansatz method, the main difficulty in the 
computation of the correlation functions lies precisely in the highly complicated structure of 
such a state and its associated excited states. Indeed, for generic, really interacting models of 
this class like the XXZ Heisenberg spin-^ chain in the regime A > — 1, the ground state is given 
as an extremely non-local action of the local fields (the spins on each lattice site here) of the 
model on a known reference state (here the ferromagnetic state with all spin up). 

This fact can be seen easily, by using the algebraic Bethe ansatz approach 0. The main 
object of this method is the so-called quantum monodromy matrix. Its matrix elements are 
operators acting in the space of states of the model, its trace (the transfer matrix) containing 
a set of conserved charges including the Hamiltonian, while the other elements lead to the 
creation and annihilation operators of the Bethe eigenstates. The quantum monodromy matrix 
is constructed as an ordered matrix product all along the chain of the so-called quantum Lax 
operators defined in terms of the local spin operators at each site of the lattice. Hence, the 
monodromy matrix itself, and the creation and annihilation operators of the Bethe eigenstates, 
are highly non-local expressions in terms of these local spin operators. In particular, let us 
mention that such operators are sums of up to 2 M terms, each being the product of M local 
spin operators all along the chain of M sites. Then, the ground state in the regime A > — 1 is 
obtained by acting N times with one of the monodromy matrix element {M = 2N for example), 
its local description being henceforth a very complicated combinatorial problem. 

Hence, although on the one hand the local spin operators satisfy a simple algebra, and on 
the other hand, the elements of the quantum monodromy matrix obey the well-known quadratic 
(Yang-Baxter) relations, the interplay between these two sets of operators is a priori quite 
involved. In particular, the commutation relations between the local spin operators and the 
quantum monodromy matrix elements do not have a simple and manageable form. This fact, 
also related to the non-local structure of the ground state (and in fact of all the other eigenstates), 
makes the computation of correlation functions very difficult. 

Indeed, until very recently, only very few models were known for which correlation function 
can be computed in an exact and explicit way. The typical examples are the Ising model, related 
to free fermions, and conformal field theories dealing with critical or massless systems in the 
continuum. 

Beyond these models, in the framework of integrable systems solvable by means of Bethe 
ansatz [||, 12 ,[l4 18 1 , related to a quantum group structure [ 19 — 22 1 and associated to an R- matrix 
solving the Yang-Baxter equation, two main approaches have been designed to deal with this 
problem. 

One of them, described in the book [12|, is based on the algebraic Bethe ansatz method |j, 12 



1£ , un applied to models defined in a finite volume (or finite lattice) . Determinant representations 



of correlation functions are obtained |j,12, 2^,^4j, containing however vacuum expectation values 
of auxiliary (quantum) "dual fields", which cannot be eliminated in general from the final result. 
Hence, explicit expressions for the correlators in the thermodynamic limit cannot be obtained 
directly from this approach. Instead, the strategy is to embed these determinant formulas in 
systems of integrable integro-difference equations from which only large distance asymptotics 
of some correlation functions can be extracted from the resolution of (matrix) operator valued 
Riemann-Hilbert problems. Note however, that for spin chains no such asymptotics have been 
determined yet, except for the free fermion points p5(|. 
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The other approach relies on the study of form factors and correlation functions of quan- 
tum integrable models directly in the infinite volume (or infinite lattice) limit. For quantum 
field theory models, it started with the study of analytic properties and bootstrap equations 
for the factorized S-matrices and form factors of integrable quantum field theories in infinite 
volume [26-28J. A typical model here is the two-dimensional Sine-Gordon relativistic quantum 



field theory. There it was realized that the set of equations satisfied by the form factors are 
closely related to the q-deformed Knizhnik-Zamolodchikov equations arising from representation 
theory of quantum affine algebras, and their q-deformed vertex operators [p9|-35|. For lattice 
models, such as the the XXZ Heisenberg spin-^ (infinite) chain, it uses the corner transfer matrix 



introduced by Baxter p7|,pq-|3q| in the context of integrable models of statistical mechanics, and 
very plausible hypothesis about the representation of the Hamiltonian of the model as a central 
element of the corresponding quantum affine algebra (here U q {sl2)) in the infinite lattice, the 



space of states being identified in terms of its highest weight modules ||, 11, 3£]. Form factors 
and correlation functions are then described in terms of q-deformed vertex operators, leading 
via bosonization to multiple integral formulas for them. 

In H we started the development of a new approach to this problem in the example of the 
XXZ Heisenberg spin-^ finite chain. In the present article we use our results |6| for the finite 
chain to obtain, in the thermodynamic limit, multiple integral representations of arbitrary n- 
points correlation functions of the XXZ Heisenberg spin-^ in a constant magnetic field. In the 
zero magnetic field case, our results agree with the one described in |8| for the massive regime, 
and with [|n]] for the massless regime. However, being grounded on the algebraic Bethe ansatz 
method, we are able to deal with the case of a constant magnetic field, a situation where the 
infinite quantum affine symmetry algebra used in |§| breaks down. 

Although developed in the context of the XXZ Heisenberg spin-^ chain, we expect that the 
general strategy we used could be applicable to other models. Our approach to the computation 
of correlation functions starts with the above general observation that the main difficulty to 
be solved is the fact that two types of operators are mixed in such quantities : the local (spin) 
operators and the (highly) non-local (in terms of the local spins) creation/annihilation operators 
of the eigenstates of the Hamiltonian. The solution to this problem is to express both types 
of operators, in a simple and manageable way, into a common algebraic structure in terms of 
which the corresponding correlation functions can be evaluated. It is this strategy that has 
been already used in fact in all cases where explicit results have been obtained for correlation 
functions : for example in the Ising model using Clifford algebra |^l|] or for the XXZ Heisenberg 
spin-^ (infinite) chain in the massive regime and in a zero magnetic field by using the U q (sl2) 
algebra and its associated g-deformed vertex operators || . Here we will solve this problem first 
for the finite XXZ Heisenberg spin-^ chain by using the Yang-Baxter algebra generated by the 
quantum monodromy 2x2 matrix, 



T(X) 



A(X) B(X) 
C(A) D(X) 



The ground state of the model can be constructed by means of algebraic Bethe ansatz by the 
successive actions of the operators B(Xk) on the ferromagnetic reference state with all spin up, 
for a particular set of spectral parameters {A^} solving the Bethe equations. 

Our method is based on the algebraic Bethe ansatz and goes along the following main steps. 
The ground state (ip g | of the XXZ Heisenberg spin-| finite chain is given as the successive 
action of the operators C(Xk) (resp. B(Xk)) (elements of the quantum monodromy matrix) on 
the ferromagnetic reference state (0| (resp. |0)), the state with all spin up, for a particular 
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set of spectral parameters solving the Bethe equations, namely, (ifi g \ = (0\Y1 C(^fc) an d 

k 

\ip g ) = Y\B(\k)\0). Our main point is that any local elementary operator Em' m can be 

k 

expressed in terms of the operators entries of the above quantum monodromy matrix || (namely 
by actually solving the quantum inverse scattering problem) as, 

, j-l A/ 

E T 3 = \\{ A ^)+D{i k ))T ejdj {i 3 ) H (A(£ k )+D(£ k j). 

k=l ' k=j+l 

Then, the elementary building blocks of the correlation functions, 

<vvi(n^ e< ) lv,ff) (L3) 

i 

can be evaluated as the scalar product of the ground state | ip g ) with the state ( ip g \ E\ l 

To compute such a quantity, we first evaluate the state (ip g l^ll^^^- There, Ml^i* 6 ^ * s 

replaced by a certain products of elements of the quantum monodromy matrix. Using their 
known (Yang-Baxter) commutation relations with the operators C(A^), it is possible to obtain 
the above wanted state as the following linear combination of states, 

(^i(n^) =5><°i n (1-4) 

^ i ' i&I keK, 

with some (computable) coefficients «j. It then remains to evaluate the scalar product of the 
ground state | ip g ) with any state in this sum. This problem was also solved in a compact form 
in (see also [^]) as the ratio of two explicit determinants. The next step is to compute the 
resulting sums in the thermodynamic limit. There, we use these determinant representations 
of the scalar products, to finally obtain multiple integral formulas, containing as integration 
measure, the determinant of the derivatives of the spectral density p characterizing the ground 
state, and solution of the Lieb equation. As a consequence, we can obtain the result also in the 
presence of a constant magnetic field. 

The article is organized as follows. After briefly describing the algebraic Bethe ansatz ap- 
proach to the XXZ Heisenberg spin-i model in section 2, the general strategy and tools of our 
method are explained in section 3. In section 4, we apply this scheme to the so-called "emptiness 
formation probability" , which is the probability to find a ferromagnetic configuration of length 
m in the massless or massive anti-ferromagnetic ground state. The computation is given first for 
a zero magnetic field in order to explain the basic techniques of the thermodynamic limit. Then 
in the next section, we derive the multiple integral formulas for an arbitrary n point correlation 
function, also for the zero magnetic field case. We generalize this result in section 6, to the 
constant magnetic field situation. 

2 The XXZ spin-| Heisenberg chain 

The Hamiltonian of the chain of finite length M [l],|2| is given by, 

M 

Hxxz = + <° y m+i + ~ !)}> ( 2 -5) 

m=l 
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and we impose periodic boundary conditions. Here a^, a = x, y, z, are the Pauli spin operators 
acting in the local quantum spin ^ space TL m at site m. The anisotropy parameter A defines the 
physical nature of the model: when A < — 1, the ground state of the Hamiltonian is ferromag- 
netic, whereas its magnetization is equal to zero when A > — 1. We shall focus our attention 
on this last domain, which itself decomposes into a massive regime (for A > 1), and a gapless 
regime (for — 1 < A < 1) in the thermodynamic limit (M — > oo). 
The i?-matrix of the XXZ model is, 



R(X,H) 



( 1 \ 

6(A,/i) c(A,Ai) 

c(A,/i) b(X,fx) 

V o o o i / 



(2.6) 



where the functions b(X, /z) and c(A, fi) are defined as 

sinh(A — /i) 



b{\n) 



c{X,n) 



sinh T) 



sinh(A — fj, + rj) ' v ' sinh(A — (j, + rj) 
The parameter r? is here related to the anisotropy parameter A of the Hamiltonian by, 

A = ^(q + q~ l ), withq = e v . 

The i?-matrix is a linear operator in the tensor product of two two-dimensional linear spaces 
V± ® V2, where each V{ is isomorphic to C 2 , and depends generically on two spectral parameters 
Ai and A2 associated to these two vector spaces. It is denoted by i?i2(Ai, A2). Such an i?-matrix 
satisfies the Yang-Baxter equation. 

Identifying one of the two vector spaces of the i?-matrix with the quantum space Ti m , one 
defines the quantum L-operator of the inhomogeneous chain at site m as, 



Lm{X,£,m) — Rom(X — 



(2.7) 



where £ m is an arbitrary inhomogeneity parameter attached to the site m. Here Ro m acts in 
Vo<g>Ti m , where Vo is an auxiliary space isomorphic to C 2 . The monodromy matrix is constructed 
as an ordered product of such L-operators: 



T(A) = Ro M (X - Cm) • • • Ro2(X - &)Roi(X - 6 



A(X) B(X) 
C(X) D(X) 



In the last formula, the monodromy matrix is represented as a 2 x 2 matrix in the auxiliary 
space Vo, whose entries A(X), B(X), C(A), and -D(A) are operators in the quantum space TL of 
the chain. 

The transfer matrix T(X) is defined as the trace A(X) + D(X) of the monodromy matrix. 
Transfer matrices commute with each other for different values of the spectral parameter A. They 
commute also with the Hamiltonian ( |2.5| ) in the homogeneous case where all £ m are equal as 
the Hamiltonian can be reconstructed in terms of the transfer matrix by means of the following 
"trace identity": 



d 

Hxxz = 2 sinh 77 — log T(A) 



+ const . 



(2i 
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Common eigenstates of the transfer matrices (and thus of the Hamiltonian ([^) in the ho- 
mogeneous case) can be constructed by successive actions of operators B(X) on the reference 
state |0), which is the ferromagnetic state with all the spins up. More precisely, the state 
-B(Ai) . . . B(\n) \ ) is a common eigenstate of the transfer matrices if the set of spectral param- 
eters {Aj}i<j</v is a solution of the Bethe equations, 

where o(A) = 1 and d(X) = \\^L\ b(\tii) are the eigenvalues of operators A(X) and D(X) respec- 
tively on the reference state | 0). The corresponding eigenvalue for the transfer matrix T(fi) is 
then, 

n n 

r(ji, {Xj}) = a(n) J] b-^Xj^) + d{n) \{ b' 1 ^, Xj). (2.10) 

3=1 3=1 

The Bethe equations can also be written in a logarithmic form: 

N 

Mp 0tai (Xj) + ^2e(Xj-X k )=2im j , l<j<N, (2.11) 
fc=i 

where rij are integers for N odd and half integers for N even. The bare momentum po to t W an< ^ 
the scattering phase 6(X) are defined as, 

M 



i . d(X) 1 v-^ ^ r l\ 



sinh(A-2 
Po(A) =zln , , 
smh(A -I- j) 

sinh(r? + A) 



k=l 

2) 



0(A) = i In 



sinh(ry — A) 



In the thermodynamic limit (M — > 00), these Bethe equations for the ground state become an 
integral equation for the quasi-particle density p in the rapidity representation (Lieb equation) 
H: 

Ptot(a) + f A K(a - p)pnt{p) dp = (2.12) 
J -A 27r 

where the new real variables a are defined in terms of general spectral parameters A differently 
in the two domains: 

a = X for - 1 < A < 1, 
a = iX for A > 1. 
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The density p is defined as the limit of the quantity M ^ a .) ; an d the functions K(a) and 
p' 0tot (a) are the derivatives with respect to a of the functions — ^ ana - Po to t(^( a )) : 

sin 2C 



Po(«) = — , -c, w ~ ( 2 - 13 ) 



Po(a) 



27T sinh(a + i£) sinh(a — i£) 
sin£ 

sinh(a + sinh(a — i|) 

sinh 2C 
27r sin(a + i() sin(a — i£) 
sinh£ 

sin(a + z|) sin(a — i|) 

M 



for — 1 < A < 1, with £ = ir/, 



for A > 1, with C = ~V, ( 2 -14) 



1 y—r £ 

with PLt( a ) = M^2Po( a - -ij)' ( 2 - 15 ) 
1=1 



where = ^ in the domain — 1 < A < 1, and (3^ = i£,k m the domain A > 1. The integration 
limit A is equal to § for A > 1, and to +oo for — 1 < A < 1. 

The solution for the Lieb equation ( [2.12 ) in the homogeneous model where all parameters 



are equal to r//2, that is the density for the ground state of the Hamiltonian (2J3) in the 
thermodynamic limit, is given by the following function [Q]: 

- 20^(f) for - 1 < A < 1, (2.16) 

P(<*) = 7T~ VT^ for A > L ( 2 - 17 ) 

K ' 2tt ^ cosh(nC) 

n=— oo v ' 

For technical convenience, we will also use in the following the solution of the inhomogeneous 
Lieb equation, that is the function, 

1 M C 
Ptot(<*) = —J2p( a -Pk-ip- (2.18) 

i=l 

It will be also convenient to consider, without any loss of generality, that the inhomogeneity 
parameters are contained in the region — £ < Im/3j < 0. 

Let us mention at last that the ground state of the XXZ model in the region A > 1 is 
degenerated in the thermodynamic limit (M — > oo), namely there are two states with the same 
energy (and characterized by the same density (2.17)), which we will call the ground state | ) 
and the quasi-ground state | ^2 } (on the finite lattice, these states possess different energy). In 
this domain, the correlation function at zero temperature is thus half of the trace on these two 
states, that is of the sum of the two corresponding matrix elements. In the domain —1 < A < 1, 
the ground state is not degenerated. 



3 Quantum inverse scattering problem and correlation func- 
tions 

In this section, we explain in more detail our general procedure to compute correlation 
functions of the XXZ chain in the algebraic Bethe ansatz framework, along the lines described 
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in the introduction. Our method is based on the study of the finite chain which has been 
performed in || . We recall here the main results of || . 

Our purpose is to compute, in the algebraic Bethe ansatz formalism, general matrix elements 
of products of local spin operators, 

( o i n (n °?) n ^ \ ° > ( 3 - 19 ) 

3 ^ i ' k 

between two Bethe states ( | \\ C((j,j) and f\ B(Xk)\ }. Here cr?% Si G {+,— ,z}, are Pauli spin 

j k 

operators at site i. The different difficulties which appear here (and which we described in the 
introduction) lead us to decompose the computation of such correlation functions into four main 
steps. 



3.1 Solution of the quantum inverse problem 

The algebraic Bethe ansatz method is based on the commutation relations given by the 
i?-matrix for the generators A, B, C, D of the Yang-Baxter algebra. The first combinatorial 
problem which appears in the study of expressions of the form fl3.19| ) is that they contain 
both local operators (spin operators) and highly non-local ones (operators B and C) and the 
relations between these two types of operators are not clear. One needs to link these two types 
of operators. One way to do it, in the spirit of the classical inverse scattering method, is to solve 
the quantum inverse scattering problem for the chain, that is to express the local quantum spin 
operators in terms of the generators A, B, C, D of the Yang-Baxter algebra. This has been 
done in ||, and we recall here the result. 

Theorem 3.1. Local spin operators at a given site i of the inhomogeneous XXX or XXZ Heisen- 
berg chain are given by, 

i-l N 

o-- = H(A + D)(£ a ) ■ B(£i) ■ J] (4 + £>)(&), (3.20) 

a=l ce=i+l 
i-l N 

<rt = H(A + D)(ta) ■ C(fc) • J] (A + D) (£«), (3.21) 

a=l a=i+l 

i-l N 

at = H(A + D)(Z a ) ■ (A -Dm ■ (A + D)(£ a ). (3.22) 

a=l a=i+l 

This reduces our problem to the mere computation of matrix elements of products of A, 
B, C, D operators in the reference state, for which commutation relations of the Yang-Baxter 
algebra can be used. 

Indeed, for any integer k and any subset {ij}x<j<k of {1, ... , iV}, with the convention 
i\ < %i--- < ik, the correlation function for spins at sites ii, ■ ■ ■ ,ik between two Bethe states 



S 



(0 I C(fii) . . . C(pb ni ) and B(X\) . . . B(X n2 ) | 0} has the following form: 

(0 | COii) • • • COinJ all °% 5 ( A i) • • • B ^ 2 ) | 0) = 

ii— 1 n\ N ri2 

= nn r1 ^ ^) • n n* -1 ^^) >< 

a=l j=l a=ik+lj=l 

ta-l 

x (0| C( M i)...C(/i ril ) • X £1 (£n) ' II + • ^ £2 (& 2 ) ••• 

o=il+l 

... J] (A + D)^) ■ X e *(£ ih ) ■ B(X 1 ) . . . B(X n2 ) |0>, (3.23) 

a=ife_i+l 

where £,-, 1 < j • < k, takes the values +, — , or z, X £j being equal respectively to C, B and 
A — D. 

3.2 Action of operators A, B, C, D on a general state 

The second step is thus to express the successive action of any product of A, B, C, D 
operators on a state constructed by action of C operators on the reference states. Action of A, 
B, C, D on such a state are well known (see for example [jl2| ), but we recall them here in a 
more convenient form for our purpose. 

The action of the operators A(X) and D(X) on the states constructed by successive actions 
of operators C(A) can be written in the following form: 

N 

n n+i 11 sinh(A fc - A a / + rj) N+ i 

(o\l[c(x k )A(x N+1 ) = Y,a(x a ,)^ Ti (0| [] C(Afc); 

fc=i «'=1 n smh(A fc - X a> ) k k ;l, (3.24) 

fc=i 

k^a 1 

N 

N N+i EI sinh(A a - X k + rj) N +i 
(0\l[C(X k )D(X N+1 ) = Y,d(X a )^ Tl (0| IJC(A fc ). (3.25) 

fc=1 a=1 nsinh(A a -A fc ) j£ 

fe=i 

The action of the operator B(X) is more complicated: 

N 

n n+i n smh (A a — A fc + 77) 

k=l 

X 



[0\HC(X k )B(X N+1 ) = Y,d(X a/ vij 

*=l «=i II sinh(A a - A fc 



k=l 



N+l 

n sinh(Aj - A a / + 77) 

N+l ,^ , j=i N+l 

^ II smh(A i -A a ^-^ (3.26) 

2=1 
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but in the case which is interesting for the computation of the correlation functions when Ajy+i = 
£k and hence d(Ajv+i) = 0, we obtain a more simple result: 



N 

N N EI sinh(A a - X k + 77) 



(0\l[C(X k )B(X N+1 ) = ]Td(A a 



, k=l 



*=i «=i n sinh ( A « - A *o 

N 

]1 sinh(Aj - A a / + rj) 

N+l 7=1 AM-1 

x E a ( A «')^ (°i n 

: '=l n S mh(A,-A a (3.27) 

3=1 

It should be mentioned that the action of B is similar to the successive action of D and A. 
Using these formulae, one can reduce expressions of the form (3.23) to sums of scalar products 
of a Bethe state with an arbitrary state constructed by successive actions of B operators on the 
reference state. 



3.3 Scalar products 



The third step of the computation is thus to find an explicit and convenient expression for 
such scalar products. Usual Bethe ansatz techniques, based only on the use of commutation 
relations, generally generate huge sums which are difficult to sum up. In ||, a direct computation 
in a new basis (.F-basis) |44j] has been performed and leads to an explicit expression for such 
scalar products as a determinant of usual functions of the model: 



Theorem 3.2. Let {Ai, . . . , Ajv} be a solution of the Bethe equations (2.9) and {^1, 
an arbitrary set of parameters. Then the scalar product, 



N N 

s N (M,{^}) = (o\ ncV;) II B ( A * 

j=\ k=l 



|0> 



,/jjv} be 



(3.28) 



can be represented as a ratio of two determinants, 

s N({^j}, {Afc}) = SW({A fc }, {Hj}) 

of the following N x N matrices T and V : 
d 



detr({/ij},{A fe }) 
detU({^},{A fc })' 



-r(// 6 ,{A fc }), 



V, 



ah 



sinh(/j b - A a ) 



1 < a,b < N, 



(3.29) 



(3.30) 



where T(fj, , {Xk}) is the eigenvalue of the transfer matrix T(Hb) corresponding to the Bethe state 
I\LiB(X k ) |0) given by $2lty. 



This result is equivalent to the scalar product formula obtained in [42|. 
When particularizing this formula in the case when the two states are equal, one obtains the 
Gaudin formula for the norm of a Bethe state: 
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(oi nc(A,) n b ^ i ° ) = sinhjV * n sinh( h A r A \\ v) ^ ^((^». 

J=i ti sinh(Ao, — \p) (331) 



where $' is a iV x iV matrix the elements of which are given by: 



^=-^ in fen^x)j- (3 - 32) 

By means of this expression for the scalar product, general correlation functions for the finite 
chain can now be expressed as sums of determinants. 

3.4 Thermodynamic limit 

The last step of our method, to obtain the general correlation functions in the infinite volume 
limit, is to take the thermodynamic limit of the expressions obtained for the finite chain. This 
has already been done in a particular case in the article jj]], where the Baxter formula for the 
spontaneous staggered magnetization in the domain A > 1 has been derived by this method. 
Here, we generalize it to any n-point correlation function. 

In the thermodynamic limit M — > oo, the Bethe equations for the ground state become the 



integral Lieb equation ( 2.12 ) for the density. In a more general way, for any C°° function / 
(7r-periodic in the domain A > 1), sums over all the values of / at the point aj, 1 < j < N, 
parameterizing the ground state, can be replaced in the thermodynamic limit by an integral 
involving the density p solution of the Lieb equation |7j: 

i N r A 

T7 E ^ = / /(«W(a) da + 0(M~°°). (3.33) 
j=l J —A 

Thus, sums over determinants will become multiple integrals. 

This properties enabled us in |7j to obtain the expression of the matrix elements of the 
Gaudin matrix (3.32) in the thermodynamic limit: 

& ab (a) = -2inM{5 ab p tot (a a ) + jjK{a a -a b )}+0(M- QO ) for - 1< A < 1, (3.34) 
$^(a) = 2TTM{5 abPtot (a a ) + -^K(a a -a b )} + 0(M- c °) for A > 1. (3.35) 

These expressions will be useful in the following to compute the determinants which appear in 
the formulae for the correlation functions in the thermodynamic limit. Finally, we will obtain 
correlation functions as multiple integrals of usual functions of the model. In the next section 
we describe in details the computation of the simplest m-point correlation function using this 
method. 



4 Emptiness formation probability 

We consider now the simplest m-point correlation function: the emptiness formation proba- 
bility, i.e. the probability to detect a ferromagnetic domain of length m in the antiferromagnetic 
ground state of the XXZ model. This probability can be expressed in the following form: 
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m 

where | ) is the ground state in the massless case and any one of two ground states constructed 
by the algebraic Bethe ansatz in the massive regime. 

Using the solution of the quantum inverse scattering problem (3.22) one can express the 
operators ^(1 — crj) in terms of the monodromy matrix elements: 

2(i-^)=n(^)+^))^) n (m)+D(t k ))- 

k=l k=j+l 

So the emptiness formation probability can be written uniquely in terms of the monodromy 
matrix elements: 

N m N 

(V\X\C{\ a )X\D(t 3 )\\B(\ a )\Q) 
r(m) = <P m ({\}) a ~ l N 3-1 N — , (4.2) 

(oinc(A a )n^(Aa)io) 

a=l a=l 

where <fi m ({\}) is the ground state eigenvalue of the corresponding product of the transfer 
matrices: 

sinh(A a - £j) 



*.(w)=nn^ 

7=1 a=l J 



Using the relation (3.25) we obtain the following action of a product of the operators D(\) 
on a state constructed by the action of the operators C(A), 

N m N+l N+2 N+m N+m 

(0\i[C(\ k )l[D(\ N+j )= Y, E ■•■ E G ol ... am (Ai...A JV+ro )<0| J] c ( x >°)> 

k=l j=l ai = l a 2 = l a m = l k=l (4.3) 

°2# a l »m^»l Om-1 fe^a 1 ,...,a m 

where the function G can be written as: 

N+j-l 

sinh(A % - - A b + 77) 

m 6=1 

G ai ...a m (Ai, . . . la N+m ) = \{ d{\ aj ) h * ai ""^ . (4.4) 

j=1 EI smh(A a3 - A 6 ) 

6=1 
6^aj ,...,aj 

To compute the emptiness formation probability one should take the parameters A a for 
a > N equal to £, a -N- It means, in particular, that the sums should be taken up to a,j = N as 
d(&) = 0. 

Now we calculate the scalar products in each term of the sum. In fact we have to calculate 
the following "normalized" product: 

N—m m N 

<oi n c(\ b ) n cfe) n £(A fe )io> 

S({A l5 . . . , \ N - m , a, • • • , {A}) = ^ „ „ • (4.5) 

(01 n c(A fc ) n 5(A fc )|o> 
fc=i fc=i 
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This quantity can be easily calculated using the representation for the scalar products ( 3.29Q 
(as one of the states in the numerator is a Bethe state) and the Gaudin formula ( 3.31 ) for the 
norm of Bethe vectors. Finally we obtain, 

s /r x \ /■ t i r x1 N TT Smh(AAr_ m+ £; ~ ^N-m+j) TT TT Smh(Afc - XN-m+j) 

S({A 1 ,...,A Af _ m ,e 1 ,..,U},{A}) = ll smm _ Cj) JUL II sinh(A fc -0) 

j>k 

AfV sinh(A a -& + i/) det^({A},{e» ( , 

M jJLi sinh ( A « " X N- m+ k + r?) det $'({A}) ' 1 ' J 

where the N x N matrix $' is the Gaudin matrix (3.32). The first N — m columns of the N x N 
matrix ty' are the same as in the Gaudin matrix but the other columns are different: 

Kb = Kb, b<N-m, 

sinhr/ 

sinh(A a - Cb+ m -N) sinh(A a - £b+m-N + v) 



Kb = ^^ : ; ■ ■ ~ ■ : : — r, b > N - m. 



As the Gaudin matrix is invertible the fraction of the two determinants in (4.6) can be repre- 
sented as one determinant: 

det*'({A},{£}) 



det $'({A}) 



det(d>'- i ({A})^({A},{e})). 



The first N — m columns of the matrix $' are the unity matrix columns: 

(V^Kab = Sab, b<N-m. 

The action of the inverse Gaudin matrix on the other columns of the matrix Vl/' can be cal- 
culated in the thermodynamic limit. Using the representation for the Gaudin matrix in the 
thermodynamic limit (3.34), (3.35) and the Lieb equation ( |2.12| ) one concludes that, 

J_ y> ~/ + sinhr? r>(M-°°) 

M thi ab P^ X b) " S inh(A a -e fc )sinh(A a -e fc + r?) + 1 h 

where p(X) is defined differently for two regimes: 

p(X) =p(X), -1 <A < 1, 
p(X) =ip(iX), A > 1. 

Thus for the corresponding matrix elements we obtain 

(S'-VU = P{K ~J" +m -\ + f) + O(M-), b>N-m. 
Mptot{Xa) 

Finally the fraction of the two determinants in (4.6) can be written in a very simple form in the 
thermodynamic limit: 

^$ { '({A}) }) = A ~Pto\^N-m+a) det S({X N ^ m+1 , . . . , A,v}, {£}) + 0(M-°° ), ^ ^ 



13 



where the m x m matrix S is: 

Sab = p(\N-m+a ~ & + |)- (4.8) 

Now, using (JO), (fi.3|), ( |4.4[ ) , (4.6) and (^), we obtain the following representation for the 
emptiness formation probability: 

N N N m 

T(m) = M m Yl sinh(g fc - 6) ^ ^ •••E^(iV,...,A am },{e})n^o 1 t(Aa J ), 

vs« M; am=lam _ i=1 ai=1 j=1 (49) 

where the function is: 

({A ai , . . . , A am }, {£}) — -p-|- . 7~7~\ ~ ; det S({X ai , . . . A am }, {£}) x 



n = 


?inh(A 


fc>Z 




771 




n 


(n. 




\fc=l 



k=j+i j (4.10) 

Here we take the sums over all the values of aj, as if two indices coincide (oj = a^, j ^ fc) the 
determinant of the matrix S vanishes. 

In the thermodynamic limit for the ground state the sums can be replaced by the integrals 
( p. 33 ) and the emptiness formation probability can be expressed as a multiple integral, 

A A 

T ( m ) = rr ■ FT / dX v ' ' I dX ™H{{\ 1: X m }, {£}), (4.11) 

11 sinh(£ fe J J 

k<l —A -A 

where A = A for -1 < A < 1 and A = -iA for A > 1. 

Thus we obtain explicit results for both regimes of the XXZ model. In both cases the 
determinant of the matrix S can be calculated explicitly. In the massless case it is the Cauchy 
determinant: 



,„ nsinhf(a-6)nsmhf(A -A 6 ) 
^=(^) ffi-Si • (4-12) 

n n sinh|(A a -a) 

a=l k=l 



2C 



The emptiness formation probability has in this case the following form: 



oo oo 



sinh|(£ fc -£z) r d\ 1 r dX m yr sinh^(A a -A 6 ) 
sinh(&-&) J ^T" J y smh(A a - \ - iQ 



k<l _ oc _ QO a>b 

mm m / j — 1 

x n n ^.rx -*a n n - & - o n - ^ 



1 lfc A = 1 1 sinhf(A a -e 



j=i yfc=i k=j+i ) (4.13) 
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In the homogeneous limit (£j = — i£/2, Vj) we obtain the following result for the emptiness 
formation probability: 



r(m) =(-1)' 



m(m + l) OO 



7 2tt' 



dA m -q sinh|(A a -A b ) 



2n LL sinh(A a - \ b - iQ 

a>b 



5=1 



sinh^-^Aj - i|) sinh m - J '(A i + i§) 



cosh™ |Aj 



(4.14) 



In the massive case the determinant of the matrix S is more complicated but can be expressed 
in terms of the Theta functions: 



det S = g n 



1 \ m j<k 
2n 



J,fc=l 



where 



n(i 



1 -Q 2n \ 2 [ n i 



,2n 



2^ -9 



2n\3 



n=l 



n=l 



m— 1 



(4.15) 



We give the proof of this formula in Appendix A. As usually in this regime we change the 
variables for more convenient ones: (3 = i£, £ = —rj. For the emptiness formation probability 
we obtain: 



k<l 



sin(/3 fe - Pi) 



ir/2 



m I m 



-tt/2 
a-l 



4/2 v =1 J a>b 



a=l \fc=l u a fc=l fc=a+l / 



^l(A a -A 6 ) 
sin(A a - X b - i() 

(4.16) 



In the homogenous limit (5j = —i(/2 we have the following result: 



n=l 



1 - g 2 " \ 2 r j 



q(m + l) -■ 



tt/2 



re=l 



.dAi 

1 27T' 



7T/2 



.d\ r 



2vr 



-tt/2 



-tt/2 



^ 2 (J? 5 n sin(A „ - a 4 - i n (4i7) 

We have shown in this section that the simplest m-point correlation function for the XXZ 
model can be expressed using the algebraic Bethe ansatz as multiple integrals of the elementary 
or elliptic functions. These results reproduce for this particular case the formulae obtained by 
Jimbo and Miwa [|, 10, Dfl. In the isotropic XXX limit (A = 1) we also obtain a formula derived 
in [^,|4(| (using the Jimbo and Miwa method). In the next section we apply the same approach 
to calculate more general correlation functions of the XXZ chain. 
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5 Correlation functions 



In this section we consider a more general case of correlation functions : the ground state 
mean value of any product of the local elementary 2x2 matrices E^ e = Si j€ iSk, € '- 

m i 

(Ainj/H) 

It should be mentioned that an arbitrary n-point correlation function can be obtained as a sum 
of such mean values. 

To calculate this product we use at first the solution of the quantum inverse scattering 
problem (3.20)-(3.22), representing the local elementary matrices in terms of the corresponding 
monodromy matrix elements: 

j-i M 

k=l ' k=j+l 

Thus we reduce the problem to the computation of the ground state mean value of an 
arbitrary ordered product of the monodromy matrix elements, 



F m{{ej, Cj}) = (f>rn ({A}) J^TJPj ' ^ ' 



where (f) m ({\}) is the ground state eigenvalue of the corresponding product of the transfer 
matrices: 



smh(A a - ij 



Now to calculate these mean values we use the commutation relations of the monodromy 
matrix elements. 

An arbitrary product of the monodromy matrix elements can be treated in a rather general 
way. At first one should consider the two following sets of indices: 

a + = {j : l<j< m, €j = 1}, card(a + ) = s', m&x jea +(j) = j^ ax , mm j£a +(j) = j' min , 
oT = {j : I <j < m, e'j = 2}, card(oT) = s, max igQ -(j) = j max , mm jea -(j) = j min . 

It should be mentioned that in a general case the intersection of these two sets is not empty and 
corresponds to the operators B(£j). 

Consider now the action of an arbitrary product on a state constructed by the action of the 
operators C(A), 

N 

(° I II c '(Afe)r ei)e ' i (A A r + i) • ..T £mje > m (\ N+m ), 
k=l 

applying one by one the formulae (3.24)-(3.26). For all the indices j from the sets a + and a~ 
one obtains a summation on the corresponding indices a'- (for j E a + , corresponding to the 
action of the operators A(X) or B(X)) or ctj (for j € a~, corresponding to the action of the 
operators D(X) or B(X)). As the product of the monodromy matrix elements is ordered these 
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summations are also ordered and the corresponding indices should be taken from the following 
sets: 

Aj ={b : 1 < b < N + m, b^ a k ,a' k , k < j}, 
A'j ={b : 1 < b < N + m, b / a' k , k < j, b ^ a k , k < j}. 

Thus the action of a product of the monodromy matrix elements can be written as the following 
sum: 

N 

(°l II C(^fc)7 1 (El)e / i (Aiv + i) . . .r £mj£ ^(AAr +m ) = G{ a . ia ' j y(\ 1 ,...,\ N+rn )(0\ C(X b ) 

k=l {aj, a 'j} 6eA m+ i (5.3) 

The summation is taken over the indices aj for j G a~ and a'j for j G a + such that: 

1 < aj < N + j, aj G A 3 -, 1 < < TV + j, a'j G A'y. 

The functions G{ a a i. }(Ai, . . . Ajv +m ) can be easily obtained from the formulae (3.24)-(3.26) 
taking into acount that A a = S,N-a for a > N: 

N+j-l 

fl sinh(A aj - X b + rf) 

6=1 

d CVi) j^g x 

J'e«- IT sinh(A % - - A 6 ) 

N+j-l 

fl sinh(A b - A a / + rj) 

6=1 J 

>< n < x ^^h ■ ( s - 4 ) 

sinh(A 6 -A a ,) 

6=1 3 

teA J+1 



Now to calculate the normalized mean value ( |5.2[ ) we apply the representation for the scalar 
product ( 3.29D and the Gaudin formula ( f3.31| ). It should be mentioned that the number of 
operators C(A) has to be equal to the number of the operators B(X), as otherwise the mean 
value is zero, and hence the total number of elements in the sets a + and a~ is s + s' = m. 
Taking into account that in (|5.2|) , for b > N, X b = (, b -N one can consider the scalar products 
appearing in the representation for the ground state mean values, 

N 

<o| n C(X b ) u B(X k )\0) 

beA m+ i k=l 

N TV ' 

<oi n c(x k ) n B(x k )\o) 

k=l k=l 

for all the permitted values of aj , a'j using the same method as for the emptiness formation 
probability. Finally we obtain: 

Fm({ej,e'j}) = sinh |g fc _^ Yl H {a j ,a' j }(^l,---,^N+m), (5-5) 
k<l {°j. a i> 
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the sum being taken on the same set of indices a,j,a'j as in (5.3). The functions H{ a . A iy({A}) 
can be obtained using (5.4) and the representations for the scalar products. It is convenient to 
introduce the following set of indices: 

{h,...,b m } = {a' j!m ,. • • , a' jLin ,a jmln , . . . , a jma J. 

One should also take into account that for the XXZ model a(A) = 1 and d(£fc) = 0. Then one 
obtains that aj < N, Vj £ a~, (otherwise the corresponding term is zero): 

^.>(w) = nS inh(A b ;-A 6i+ ,) n n**^-&-*) n u 
x n ( n sinh ( A % n ^ma^ - &) ) ( det m ({^» + o(^~°°)) , 

iea-\fc=i fe=i+i / V (3.6) 

where the mxm matrix M({&&}) is slightly different in comparison to the case of the emptiness 
formation probability, 

k > N, Mik = -6b,-N,k, 

f>{\ ~ £fc + \ ) 



k < N, Mi 



Ik 



Ptot (Aft 



The sum in ( |5.5| ) can be rewritten in a more simple way if one takes into account that the 
function H{ a j0 '.}({A}) defined by (5.6) is equal to zero if bj = b^, j ^ k (as the determinant 
vanishes in this case), or if a'- > N + j: 



N+m N+m N N 



11 v« w 6i=1 6s , =1&s , +i=1 6m=1 

The sum over 1 < < iV is just a sum over the rapidities in the ground state and can be 
replaced by integrals as in the case of the emptiness formation probability: 



N ^ 

Y,f(X a )= j d\p tot (\)f(\) + 0(M-°°). 

0=1 r 



The contributions of the terms with a'- > N can be rewritten as integrals over the contours 
r a /,_jv surrounding the pole of the corresponding density function p{\ — C a '^-N + S) m * ne P°i n t 
A = C a '.-N- The residues of the density function p{\ — £ + 3) (for both regimes) in this points 



3 

are: 



2™ Res^(A-e+|)) 



= -1. 

The other points should be outside the contour r o /_jy. The matrix M({6/%}) then should be 
replaced by the same matrix SdA}, {£}) as in the representation for the emptiness formation 
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probability. Finally for the correlation function one obtains: 

k<l \-A J / V-A J / 

A A 

J d\ s , +v ■ ■ J d\ m H{ €j ,e'.}(\i, ...,\ m ) + 0(M- 
-A -A 

where the function H^ € ^ e /y(Xi, . . . , A m ) is defined as, 

(— l) s ' i _1 m 

%,, ; ,({A}) = nsmh ^ fc J _ A;+r/) n n ^ n ^ -^\- 



(5.8) 



j£a~ \k=l k=j+l 



(j-1 m \ 

[] sinh(^- - 6 - r?) sinh(^ - &) det S({A}, {£}), (5.9) 
k=i k=j+i J 

where 

Sik = p{k ~ik+ 
and the parameters of integration are ordered in the following way: 

{Al, . . . A m } = {fJ-j' , • • • , tfjl . , ^j min , • • • , /^ max }- 

J max J mm 

Consider now separately the two regimes of the XXZ model. In the massless regime 77 = — i( 
is imaginary, the ground state rapidities A are real and the limit of integration is infinity A = 00. 
In this case we consider the inhomogeneity parameters £j such that > Im(£j) > — £. The 
function H^ e . e i }(Ai, . . . , A m ) for all the arguments Aj in the region > Im(Aj) > — £ has only 
simple poles in the points Aj = Hence the sums of integrals in (5.8) can be rewritten as one 
integral on a displaced contour: 

co— iC, 



00— i( 

Finally for the correlation functions in the thermodynamic limit one obtains the following 
result in this regime: 



00— if 00— i( 



00 



Fm({tj, = rr ■ vJl FT / ^ 1 " ' / I ^K'+v ■ ■ / d\ m H^ e . e ry(\i, ... , A m ). 

11 smn lS/c Q) J J J J 3 (x,m\ 

k<l -00-iC -00-iC _ °° _ °° v u - ±u ; 

Now one can rewrite this result using the corresponding representations for the determinants 
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of the matrix 5({A}, {£}): 



mm I j—1 m 

x iin sinh ^ A n ri sinh (^-^-<) n ^^-^) 

a=lfc=l C v a jea~ \k=l k=j+l 

(j—1 rn \ 

n sinh^ - & + »o n sinh (^ - &) • ( s - n ) 
fc=l k=j+l j 

In the homogeneous limit (£j = —i(/2, the correlation function i^ m ({ej, e' }) has the 
following form: 



7T\ 2 rr f dXj -n f d\j t~T smn £ l A a _ A fc) 



=(-!)■ (-f) 11 / : ii / : ii 



. . .. - , . ., 2vr 11 sinh(A a - A 6 - <) 



X 



-p-j- sinh J_1 (^ - i|) sinh m_J (/ij + i|) 
,„,,„/-'( ,/ j- { v±-\ "'-■'/ ,,' _i_ 

x 



cosh m |/^- 



sintf-^/i'. +34)sinh m ^(/i' / s 

n — — — ■ (5 - i2) 

j<=a+ i J 



These results agree exactly with the ones obtained by Jimbo and Miwa in [10], taking into 
account that their Hamiltonian is obtained from our one by the transformation UH&U^ 1 = 



M 



-H_ A ,U = Uf =1 ^- 

Similarly in the massive regime the parameter r/ = — £ is real, the solutions of the Bethe 
equations corresponding to the two ground states are imaginary A a = — ia a , and A = —in/2. 
The inhomogeneity parameters are chosen in such a way that > Im(/5j) > — £, £ = —if}. 
Taking into account that the function H^ e ,}(Ai, . . . , A m ) is 7r-periodic for all the arguments 
a a and that in the region > Im(a a ) > — £, — ir/2 < Re(a a ) < 7r/2, it has only simple poles 
in the points a a = (5k one can rewrite the sums of integrals in the same way as in the massless 
situation: 

da y 

■w/2-iC, 




20 



Finally using the formula for the determinant of the matrix S ( f4.15| ) in this case one obtains: 



F m({ej,e'j}) =9mYl 



k<l 



sin(/3 fc - fii) fJL 



, TT/2-iC 



■k/2 



-7r/2-iC 

0-i 



i=s'+i. 



^i(A a -A 6 ) 



•tt/2 



27r sin(A a - Aft - <) 



a>6 



x iifi xa n (n^-A-io n 

a=lfc=l u a HK > j(z a - \k=X k=j+l 

(j—X m \ I m 

H sin(^- - & + iC) n sin(/4 - tf 2 £>j - # 
fc=i fc=i+i / y=i 

where {Ai, . . . A m } = {/4 , ■ • • , . , /% min , • • • , ^j max } and 




(5.13) 



min 

1 



9m = 



n=X 



+ 



-,2n 



n=X 



m—X 



This result is equivalent to the formula [||,[Tl[] obtained as a solution of the q-KZ equations. 

In the homogenous limit (5j = —iC/2 we obtain the following result for the correlation 
function: 



1 



^({^•})=n(i+ 
n 



2« N 2r 



-,2n 



2n\3 



n=X ' n=X 

>6 sin(A a -A b -<) . e „ 



, ir/2-iC tt/2 
x(m+i) , s' ' , j, m L 

/ £ n 



i=i 



7r/2-iC 



i=s'+l_ 



•tt/2 



-j-j- sin-' — i|) sin m 3 '(/ i i + *2> 



n 



sin?- V 7 - + 34) sin m ~Vi + »§) 



tf 2 + 



(5.14) 



Thus for both regimes the results for the correlation functions obtained by Jimbo, Miwa and 
collaborators 11] as solutions of the q-KZ equation can be reproduced using the algebraic 

Bethe ansatz approach. We have shown also that the main difference between the two regimes 
is the determinant of the matrix the elements of which are the density function taken at the 
corresponding values of the spectral parameter. In both cases this determinant can be calculated 
in terms of elementary or elliptic functions, which leads to the different representations for the 
two regimes. 

It should be mentioned that to calculate the correlation functions one can also use the action 
of the monodromy matrix elements on a state constructed by the action of the operators B(X). 
It leads to a similar but different representation for the correlation functions. 



6 External magnetic field 

In this section we show how the previous results can be generalized for the XXZ spin chain 
in a constant external magnetic field: 



H, =H 



XXZ 



hS z , 



(6.1) 
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where S z is the third component of the total spin, 



M 
3=1 

The third component of the total spin commutes with the Hamiltonian of the XXZ model and 
also with the transfer matrix and so the eigenstates remain the same as in the case of zero 
magnetic field. However the ground state changes and the corresponding Lieb equation changes 
also. In the massless regime the ground state changes for any value of the magnetic field. The 
Lieb equation now has the following form, 

Ah 

Ph (X)+ I dpK(X-p)p h ( f i) = ^p' (X), (6.2) 

-Ah 

where the Fermi momentum is defined by the following integral equation for the excitation 
energy: 

A,, 

e h (X)+ J dnK{\- n)e h (ji) = e (X), e h (A h ) = 0, (6.3) 

-Ah 

where the bare energy £o(A) is, 

e = h - 2sinCPo(A). 

In the massive regime the ground state changes only for the magnetic field greater than the 
critical value equal to the gap width. Indeed the equations defining the Fermi momentum, 

A,, 

e h (X)+ [ dpK(X- p)e h (p) = e (A), e h (A h ) = 0, (6.4) 



-A h 

with the bare energy 

Eq = h — 2sinh£po(A), 
have no solution for < tt/2 if the magnetic field h < h c , where the critical field is: 

/^smhc f; 

7T cosn nC 

n=—oo 

So if the magnetic field is under its critical value, the ground states, the density function p(A) 
and the zero temperature correlation functions do not change. If the magnetic field is greater 
than h c there is a solution for the Fermi momentum A^ and the density function is given by 
the equation Q6.2j) with the kernel K(X) and the function p' (X) corresponding to the massive 
regime. It should be mentioned that in this case there is no more gap in the spectrum, and the 
XXZ model with A > 1 and h > h c is massless. 

In general these equations cannot be solved explicitly as in the case of zero magnetic field, 
but as the integral operators are rather simple the solution can be obtained numerically with 
any given precision. 
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The function Ph(X) has only one simple pole in the region > Im(A) > — £ in the point 
A = — i£/2 and its residue is, 



2ni Res I ph(X) 



1. 



The correlation functions can be calculated using exactly the same method as in the previous 
section. It should be mentioned that in this case only the homogeneous limit is interesting. Using 
the same arguments one can show that in the regime |A| < 1 the general correlation function 
can be written as follows, 



F m (/ l ,{6„ e ;.})=n J + J +J \dXj n J dA m 



(-1)*' 



.n(- 



sinb 3 (/Uj - i-) sinh m + . Q 



A h k>l 



sinh(A fc - A; - i() 



H (smb*- 1 ^ + 3i|) sinli'" ^ + det^({A}), (6.5) 



where {Ai, . . . X m } = {//•, , . . . , /// , fJ>j miu , • • • , Mimax} an d the matrix elements m x m matrix 

J max J lYiin 

'S'hafe = Ph,b{^a), 

are defined as solutions of the following integral equations: 

A,, 

Ph,bW + J dpK(X - fj,)p h>b (fi) 

-A h 



d b 



2Tr(b-l)\dX b 



Po(A). 



(6.6) 



In the case A > 1, h > h c the result has the following form, 



/ -7T/2 TT/2-iC A,A 



F m (/ l ,{e J , e ;.})=n 



+ 



+ 



A,, 



\-Aft -7r/2-iC tt/2 / 



(iAo 



(-i)' 



sin^' ^ -A)sin m X 



, ; •' II sin(A fc - Xi - i() 



sin^ 1 ^ + 3tJ) sin m ^(^ + J) ) detS„({A}). 



-C 



(6.7) 



The matrix elements of the matrix are given also by the integral equations ([O]) with the 
kernel K{X — p) and bare momentum p$(X) corresponding to the regime A > 1. 

The determinants cannot be calculated explicitly like in the case of zero magnetic field but 
however the density functions and hence the determinants can be calculated numerically from 
the corresponding integral equations. 



Acknowlegement. We would like to thank A. Izergin and N. Slavnov for useful discus- 
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Appendix A 



In this appendix we compute the determinant of the matrix (Sij)i<i ; j< m , with Sij = p{Xi — 
j3j — of the density in the massive case. In the domain A > 1, p is an elliptic function with 
periods ir and 2iQ (q = e - ^), which can be written in terms of Theta- functions: 

j +°° e 2i\n 2 q 2n \2 fl 3 (\ q\ 

pW = 2^ cosh< = 2^^\l + q 2n ) M\q)' ^ 

n=— oo n=l / 

To compute the determinant det m 5, let us first consider it as a function / of the variable Ai. 
/ is thus an elliptic function of Ai of order 2m and with the same periods as p. An irreducible set 
of poles is . . . , P m , Pi +iC) • • • > Pm + *C}> an d A2, • • • , A m are obviously zeros of /. Note that, 
for any x, f(x + i() = —f(x), hence A2 + i(, . . . , A m + i( are also zeros of /. Up to congruence, 
there remains also only two other zeros which differ by i(, say xq and xq + iQ. Since the sum of 
zeros of an elliptic function is congruent to the sum of its poles, it follows that x$ is congruent 

mm mm 

either to ^ Pi — ^ Aj or to Yl Pi ~ Ai ~~ § • Actually it is congruent to the second expression, 

i=l i=2 i=l i=2 

as it will be shown latter. 

Let us now consider the function 



n^MAi-A,) 



n^i(Ai-^- 

5 is an elliptic function with periods ir and 2i£, of order 2m, which has the same poles and zeros 
as /. By Liouville's theorem, f /g is a constant. Note at this stage that if xq was taken to be 

m m 

congruent to Pi ~ Ai, a similar argument would have lead to the fact that / should be 

i=l i=2 

equal, up to a multiplicative constant, to, 

m 

n^i(Ai-A,) 

£ ^E(Vft). 

n^i(Ai-^) s=i y 

which is obviously not true because the periods do not coincide. 

The same procedure for the variables A2 , • • • , A m , and similarly p\ , . . . , (3 m , leads to the 
following formula for det m S: 

UMXj-x^MPk-Pj) m 

det m ~S = Cj— ^(^(A,- -Pj)), (A.2) 

n M^j-pk) s=i 7 
j,k=i 

where C m is a constant which does not depend on Aj, Pj, 1 < i, j < m. 

A recursion relation for C m can be obtained by taking the residue of the two members 
of ( |A.2j ) at the pole A m = P m : 

I 00 



C m = i—T\{l-q 2n ? C n 



^m— 1 ■ 

7T " " 

71=1 
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The determination of C\ is straitforward 



oo 1 



l-0 2n \2 



n=l 

and thus we obtain for the constant C m : 



+ 2n 



Cm. — I 



I ^ m T r ^ ]_ — ^ 2 |- „ , , r _ „-| m—1 



V2vr 



„2n\3 

^ - r 

n=l ' ^ n=l 
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